# Standard data
df <- readRDS("./data/dataset/df.RDS")
tr <- read.tree("./data/tree/tree.treefile")
pclustering <- readRDS("./data/asr_clustering/blbli_asr_clustering_df.RDS")
df <- df %>% left_join(.,pclustering)
# Matrices
## KPC Plasmids
KPC_info_clusters <- readRDS("./data/KPC_plasmid/KPC_containing_clusters_mat.RDS") %>% as.data.frame
KPC_containing_clusters <- KPC_info_clusters %>% `colnames<-`(paste0(colnames(.),"_KPC")) %>% {ifelse(.=="KPC Plasmid",1,0)} %>% as.data.frame%>% select_if(colSums(.)>0)
KPC_clusters <- colnames(KPC_containing_clusters)
Non_KPC_clusters <- KPC_info_clusters %>% `colnames<-`(paste0(colnames(.),"_non_KPC")) %>% {ifelse(.=="Non-KPC plasmid",1,0)} %>% as.data.frame %>% select_if(colSums(.)>0)
df <- left_join(df,KPC_containing_clusters %>% mutate(isolate_no = rownames(.))) %>% left_join(Non_KPC_clusters %>% mutate(isolate_no = rownames(.)))
Gain & Loss Figure
traits_of_interest <- c("blbli_dich_num",KPC_clusters)
trait_parent_child <- pbmcapply::pbmclapply(traits_of_interest,FUN=function(x){
asr(df,tr,tip_name_var = 'isolate_no',pheno=x,model='MF') %>% .$parent_child_df
},mc.cores=5) %>% `names<-`(traits_of_interest)
synchronous_changes_w_AA552 <- pbmcapply::pbmclapply(traits_of_interest %>% subset(.!="AA552_KPC"),FUN=function(x){
phyloAMR::synchronous_detection(comparitor_parent_child_df = trait_parent_child[[x]],trait_parent_child_df = trait_parent_child[['AA552_KPC']],node_states = 'joint')
},mc.cores=5) %>% do.call(rbind,.) %>% mutate(comparitor = traits_of_interest %>% subset(.!="AA552_KPC"))
synchronous_changes_w_AA552 %>% subset(synchronous_transitions_num>0)
AA552_gain <- trait_parent_child$AA552_KPC %>% subset(gain ==1) %>% .$child
AA018_loss <- trait_parent_child$AA018_KPC %>% subset(loss ==1) %>% .$child
AA274_loss <- trait_parent_child$AA274_KPC %>% subset(loss ==1) %>% .$child
st258 <- ggtree(tr)
tree_obj <- st258 +
geom_nodepoint(aes(fill = "AA552 KPC plasmid gain",shape="AA552 KPC plasmid gain",subset=(node %in% AA552_gain),size="AA552 KPC plasmid gain",legend=F)) +
geom_tippoint(aes(fill = "AA552 KPC plasmid gain",shape="AA552 KPC plasmid gain",subset=(node %in% AA552_gain),size="AA552 KPC plasmid gain",legend=F)) +
geom_nodepoint(aes(fill = "AA274 KPC plasmid loss",shape="AA274 KPC plasmid loss",subset=(node %in% AA274_loss),size="AA274 KPC plasmid loss",legend=F)) +
geom_tippoint(aes(fill = "AA274 KPC plasmid loss",shape="AA274 KPC plasmid loss",subset=(node %in% AA274_loss),size="AA274 KPC plasmid loss",legend=F)) +
geom_nodepoint(aes(fill = "AA018 KPC plasmid loss",shape="AA018 KPC plasmid loss",subset=(node %in% AA018_loss),size="AA018 KPC plasmid loss",legend=F)) +
geom_tippoint(aes(fill = "AA018 KPC plasmid loss",shape="AA018 KPC plasmid loss",subset=(node %in% AA018_loss),size="AA018 KPC plasmid loss",legend=F)) + scale_fill_manual(name="Gain/Loss Events",values = c("AA552 KPC plasmid gain" ="#91796A","AA018 KPC plasmid loss"="#97C466","AA274 KPC plasmid loss"="blue"),guide = guide_legend(title.position = "top",label.position = "bottom",ncol=1,order=1)) +
scale_shape_manual(name="Gain/Loss Events",values=c("AA552 KPC plasmid gain"=22,"AA018 KPC plasmid loss"=24,"AA274 KPC plasmid loss"=23),guide = guide_legend(title.position = "top",label.position = "bottom",ncol=1,order=1)) + scale_size_manual(name="Gain/Loss Events",values = c("AA552 KPC plasmid gain" =20,"AA018 KPC plasmid loss"= 16,"AA274 KPC plasmid loss"=16),guide = guide_legend(title.position = "top",label.position = "bottom",ncol=1,order=1))
A. KPC Plasmid content
